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Abstract 



We present a phase-field model for the dynamics of the interface between 
two inmiscible fluids with arbitrary viscosity contrast in a rectangular Hele- 
Shaw cell. With asymptotic matching techniques we check the model to yield 
the right Hele-Shaw equations in the sharp-interface limit and compute the 
corrections to these equations to first order in the interface thickness. We 
also compute the effect of such corrections on the linear dispersion relation 
of the planar interface. We discuss in detail the conditions on the interface 
thickness to control the accuracy and convergence of the phase-field model to 
the limiting Hele-Shaw dynamics. In particular, the convergence appears to 
be slower for high viscosity contrasts. 

Copyright 1999 The American Physical Society 
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I. INTRODUCTION 



The dynamics of morphologically unstable interfaces is a major problem in nonequilib- 
rium physics from both fundamental and applied points of view. Relevant examples of those 
are dendritic growth, directional solidification, flows in porous media, flame propagation, 
electrodeposition or bacterial growth M. The so-called Saffman- Taylor problem has played 
a central role in this context because of its relative simplicity both experimentally and in 
its theoretical formulation |IjJ2]. It deals with the motion of the interface between two in- 
miscible fluids within a Hele-Shaw cell. Due to the highly nonlinear and nonlocal nature of 
the interfacial dynamics of such systems, analytical understanding is scarce and restricted 
to high viscosity contrast 0, so in general one relies mostly on numerical work 

From a mathematical point of view, such systems are referred to as moving boundary 
problems. In practice this implies that one has to keep track of the interface where boundary 
conditions are applied, and solve a (linear) problem in the bulk which determines in turn the 
motion of the boundary. This kind of problem has traditionally been addressed in terms of 
boundary integral methods which reduce the dynamics of the interface to integrodifferential 
equations. The numerical integration of these equations is quite involved though, particu- 
larly for long times, due to stiffness and numerical instability of the equations. In the case of 
Hele-Shaw flows, boundary integral methods have succesfully been applied although 
quite sophisticated algorithms have usually been necessary |§. 

Recently, the so-called phase-field equations have been proposed in the context of solidifi- 
cation problems as a different approach to the interface dynamics []T~Q| p~9f| . In the spirit of the 
well known time-dependent Ginzburg-Landau models |2(J|| , the method avoids the tracking 
of the interface by introducing an auxiliary field (analogous to an order parameter) which 
locates the interface and whose dynamics is coupled to the other physical fields through an 
appropriate set of partial differential equations. In this way, there is no boundary condition 
to explicitely apply at the interface and the whole system is treated as bulk. 

This method introduces a mesoscale e which is not present in the original macroscopic 
equations and gives a finite thickness to the interface. The equations are then chosen in 
such a way that the original bulk equations and boundary conditions are recovered in the 
e — > limit. Therefore the phase- field equations for a given model are not intended to 
describe the true mesoscale physics of the system, and are then not unique. In fact there is 
considerable freedom in choosing a particular form of them, with criteria of either numerical 
efficiency and convergence [13[] or other physical criteria such as thermodynamic consistency 
|I4| . In any case, the nature of the phase-field approach is completely different from the 
sharp-interface models and therefore the actual numerical advantages and limitations of both 
are also quite distinct. This makes the two approaches complementary and competitive in 
different physical situations. A remarkable advantage of the phase-field approach is that it 
is much simpler to implement satisfactorily from a numerical point of view. On the other 
hand, the phase-field approach is usually more amenable to generalization, in the sense that 
it allows to introduce variations and new elements without any major modification of the 
numerical scheme, for instance in the treatment of fluctuations, liquid crystals and other 
complex fluids 0. Finally, the phase-field approach can handle very naturally situations 
where the sharp interface model is not appropriate, such as for instance topology changes 
like interface pinching leading to the breakup of bubbles. 
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In this paper we introduce a phase-field model for Hele-Shaw flows with arbitrary vis- 
cosity contrast (or Atwood ratio) c = . Although in the high contrast limit c = 1 the 
Hele-Shaw dynamics is quite analogous to the one-sided solidification problem (in the ap- 
propiate approximations ||), the arbitrary viscosity contrast case has been shown to exhibit 
quite different dynamics than solidification problems, and has in fact opened some inter- 
esting questions, particularly concerning the sensitivity of finger competition to viscosity 
contrast |^-[7|J2T|] and the long time asymptotics of the low viscosity contrast limit 0. 

The model presented here is inspired in the vortex-sheet formulation of the problem 

in which the relevant dynamic variable in the bulk is the stream function. Similar ideas 
have previously been applied to describe physically diffuse interfaces in the context of steady 
state selection in thermal plumes P^] . Usually phase-field models are naturally suited to 
symmetric situations (two-sided models). The present case of Hele-Shaw flow is no exception 
and becomes most efficient for c = 0. Finite c can also be handled but the model becomes 
computationally inefficient in the limit c — > 1, since this limit must be taken formally after 
the e — > limit. A phase-field model for this one-sided case must differ essentially from the 
one presented here, such as in the spirit of Ref. fll9| . 

The layout of the rest of the paper is as follows: in Sec. |II A| we recall the Hele-Shaw 
macroscopic equations in terms of the stream function, whereas in Sec. |lllj| we present our 
phase-field equations. We then show in Sec. [HJ that the phase-field equations reduce to the 
macroscopic ones in the sharp- interface limit. Deviations from that limiting behavior are 
derived from the phase-field equations themselves to first order in the interface thickness 
in Sec. [TV], and their effect on the linear regime is computed in Sec. [V]. Finally, a brief 
summary is given in Sec. [VT]. 



II. THE MODEL 



We consider the general case of an interface with surface tension a between two fluids 
with distinct viscosities (/ii, /z 2 ) and densities (pi, p 2 ) moving in a rectangular Hele-Shaw cell 
of width W (x-direction) and gap b (^-direction), under an effective gravity g e ff (negative 
y-direction) and with an injection velocity V^, (positive y-direction). Label 1 (2) corresponds 
to the upper (lower) fluid. 



A. Macroscopic equations 



Darcy's law is assumed to hold for each fluid, thus defining a certain velocity potential 
in each bulk, but not on the interface. In contrast, the bulk incompressibility and the 
continuity of normal velocities on the interface allow us to define its harmonic conjugate, 
the stream function ip, even on the interface through u x = d y ip, u y = —d x ip, where u x , u y 
are the x, y components of the fluid velocity field u. Then Darcy's law results in a Laplace 
equation for the stream function (potential flow) and a certain jump for the tangential fluid 
velocities on the interface, whose value takes into account the Gibbs-Thomson relationship. 
The fact that the stream function is continuous at the interface makes the use of this variable 
particularly convenient. The Hele-Shaw equations in stream function formulation |3J] can be 
written in dimensionless form as: 
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V 2 ip = 



(2.1) 



Vv(0 + ) - Vr(O-) = -7 - <# r (0 + ) + Vr(O-)] (2.2) 

^(0 + )=^ s (0-) = - Wn , (2.3) 

which constitute the Hele-Shaw equations. Here r is a coordinate normal to the interface 
and with origin on it, positive in fluid 1 (0 then means on the interface coming from each 
side), s is arclength along the interface and such that the unit vectors satisfy s x f = x x y, 
the subscripts stand for partial derivatives except for v n (s), which is the normal velocity of 
the interface, positive towards fluid 1, and 

M = Bk s + y ■ s, (2.4) 

with k,(s) the interface curvature, positive for a bump into fluid 2. The dynamics are 
controlled by the two dimensionless parameters 



B = *1 c ^ 1-/i2 f 2 5) 

12W 2 [V 0O (ix 1 -ix 2 )+g eff ^(p l -p 2 )\ ' /X1 + M2' 



We will not be interested in negative values of B (stable configuration) nor c (mirror im- 
age interface of — c). So B is a dimensionless surface tension, and can be understood 
as the ratio between the capillary (stabilizing) force and the driving (destabilizing) force 
(injection+gravity), and c is the viscosity contrast, which is so far completely arbitrary: 
< c < 1. This corresponds to having set ourselves in the frame moving with the fluid 
at infinity (or, equivalently, with the mean interface) and taken W as unit length and 

u* = cVoc + ^//SferSi as unit velocit y ( see ID- 
Note that Eqs. ( |2.1| , |2.2| ) can be written together as 

V 2 ^ = -w, w = { 7 (s) + c# r (0+) + Vr(0")]}<S(r) (2.6) 

where 5(r) is the Dirac delta distribution and w = z ■ (V x u) is the fluid vorticity, which is 
confined to the interface. 

B. Phase- field equations 

We put forward the following phase-field model for the above equations with 6 being the 
phase field: 

e ^ = v2 ^ + c ^ ■ ( ^ ) + \^r {m _ ° 2) (2 - 7) 

e 2 ^- = f(0) + e 2 V 2 9 + e 2 n(9)\V9\ + e 2 z ■ (V^ x W) (2.8) 
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where f(9) = 9(1 - 9 2 ), and = s(9) ■ {BVk{9) + y), k(9) = -V • f{9), with f{9) = ^ 
and s(#) = r(9) x 2, together with the boundary condition 

0(3/ - ±00) = ±1, (2.9) 

so that = +1(— 1) corresponds to fluid 1 (2). 7(0), k(0) are functionals which generalize 
the magnitudes defined above for the interface, now to any level-set of the phase-field. 
If we leave the two last terms aside, Eq. Q2.8Q is the Cahn-Hilliard equation for a non- 



conserved order parameter or model A (without noise) in the classification of Ref. [ 20 1 of 
time-dependent Ginzburg-Landau models. The field in this model is known to relax towards 
a kink solution of a certain width in a short time scale, and then to evolve to minimize the 
length of the effective interface according to Allen-Calm law (i.e. with normal velocity 
proportional to the local curvature). The factor multiplying the laplacian has been choosen 
to be e 2 for the kink width to be 0(e), so that e can be considered the interface thickness, 
i.e., the small parameter in the asymptotic analysis that will be performed in next section. 
On the other hand, the e 2 factor in the time derivative ensures that the relaxation towards 
the kink is much faster than the evolution of the interface. Notice that model A describes the 
relaxational dynamics of a non-conserved order parameter, whereas our problem is actually 
non-relaxational and strictly conserved (mass consevation and inmiscibility) . The other two 
term in the phase-field equation will correct this apparent contradiction. In order to cancel 
out the local Allen-Cahn dynamics of the interface which is buit in model A, we add the 
term e 2 /t(0)|V0|. It will be shown that such term cancels out Allen-Cahn law by giving rise, 
to leading order, to an identical contribution but with opposite sign. With these elements 
so far, our phase-field relaxes to a kink profile located along an arbitrary interface which 
(if sufficiently smooth) remains almost completely stationary, regardless of its shape. This 
is because the dynamical effect of surface tension associated to the Ginzburg-Landau free 
energy has been removed (up to first order) and the interface has not yet been coupled to the 
fluid flow, represented by the stream function. This coupling is achieved by adding the last 
term in Eq. Q2.8| ), which stands for —e 2 u- V0 and thus sets the phase-field — and therefore 
the interface — in the frame moving with the fluid velocity u. This term restores the fully 
nonlocal dynamics of the Hele-Shaw model. In particular it yields the continuity of normal 
velocities Eq. (|2.3|) and reintroduces surface tension, which is contained in the dynamical 
equation for the stream function through 7(0). 

As for Eq. (|2.7|) , its right hand side is intended to reproduce Eq. (|2.6|) , and therefore 
also Eqs. ( |2.1[ ) and ( |2.2p . If the phase- field has a kink shape, 1 — 9 2 is a peaked function 
which, when divided by e, gives rise to the delta distribution for the vorticity. However, this 
only accounts for the 7 in the weight of the delta. The part proportional to the viscosity 
contrast c must be put apart as the cV • (0V0) term because of the non-local character of 
Vv(0 + ) +-?/v(0~). Finally, the time derivative is multiplied by e to recover the laplacian (and 
not diffusive) behavior of the Hele-Shaw flow in the sharp- interface limit. 



In spite of important differences, the proposed phase field equations Eqs. (|2.7| , |2.8|) contain 
certain similarities to the problem of a thermal plume in a Hele-Shaw cell under gravity P2| . 
In such a problem there is only one fluid heated from the centre of the channel. The heat 
diffuses towards the lateral walls, but the temperature profile is not linear, since the fluid 
density and viscosity decrease with temperature, so that the fluid in the middle of the 
channel raises because of buoyancy. As a result a so-called plume of hot fluid with a shape 
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similar to the Saffman-Taylor finger, with a stationary upwards velocity and a width close 
to 1/2 is formed. Outside the plume the fluid is colder, and the transition between the two 
zones is relatively abrupt, so that one can think in terms of an interface of a certain small 
thickness. Thus, the equation for the phase field Eq. Q2.8|) could be thought as a diffusion 
equation for the temperature in a thermal plume. However, the available equations for that 



problem hold only for the steady state [22], whereas our phase-field model is intended to 
describe the whole dynamics. Generalization of the thermal plume equations to include 
the dynamics is not trivial for non- vanishing viscosity contrast. As a matter of fact, Ref. 
2"2| must restrict itself to low viscosity contrasts — as it is the case in thermal plumes — , 
whereas we formulate the model for arbitrary viscosity contrast. An interesting difference is 
the term e 2 K,(6)\V6\ cancelling out Allen-Cahn law. The absence of that term in the thermal 
plume equations does not prevent the Hele-Shaw steady state equations to be recovered in 
the sharp-interface limit because of the lower power of e used in the u ■ V6 term, but then 
Allen-Cahn law arises in the corrections at first order in the interface thickness. In contrast, 
by means of this e 2 K(0)|V#| term we achieve cancellation of the Allen-Cahn law even in such 
corrections, as we will see in section [TV]. Finally, another major difference in the case of 
thermal plumes is the absence of surface tension. 

III. SHARP-INTERFACE LIMIT 

In order to analyze the small-e behavior of the phase- field equations, Eqs. ( p.7||2.8|) , we 
expand their fields in powers of e. The expected abrupt variations of these fields through the 
interface will make it necessary to perform two different expansions. In the interface region 
(inner region) we rescale the differential operators appearing in these phase-field equations 
by rewritting them in terms of the streched normal coordinate p = r/e (see Appendix). The 
expansions in the inner region will be matched order by order in powers of e to those in the 
outer region (in the bulk far from the interface) where the coordinates are not rescaled. The 
outer and inner expansion are written respectively as 

a(r, s,t) = a (r,s,t) + ea 1 (r,s 1 t) + e 2 a 2 (r,s 1 t) + ... (3.1) 

A{p, s, t) = A {p, s, t) + eA x {p, s, t) + e 2 A 2 {p, s, t) + ... (3.2) 

where capital letters denote fields written in terms of the rescaled coordinate. This results 
in the following matching conditions: 

A (p,s,t) = a (0 ± ,s,t) 

Ai(p, s, t) = a 1 (0 ± , s, t) + pa 0ir (0 ± , s, t) as p — > ±oo (3.3) 

A 2 (p, s, t) = a 2 (0 ± , s, t) + pai jr (0 ± , s, t) + ya ,rr(0 ± , s, t) 



And therefore: 
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A 0jP (±oo, s, t) = A ljPP (±oo, s, t) = ... = 0, 
Al, p (±oo, s, i) = a 0ir (0 ± , s, t) 

A-2, P {p, s,t) = a l r (0 ± , s,t) + pa 0irT (0 ± , s,t) as p — > ±00 



(3.4) 



In practice, one does not find explicit solutions for the fields, but some set of equations 
for them. A sharp-interface model for the small-e dynamics of the phase-field equations, 
Eqs. (127^ , F8|) , is then given by the set of equations obeyed by the outer fields: Those 
obtained at lowest order in the interface thickness e (O(e )) constitute the e — > limit of 
the phase-field model, which we carry out in this section; whereas those obtained up to 0(e) 
represent what we will call (following Karma and Rappel |13[) a 'thin- interface' model, a 



model keeping finite interface thickness effects, such as the one derived in Sec. IV 



A. Outer equations 

Straightforward substitution of the outer expansion Eq. ( |3.1| ) in the outer equations, 
Eqs. ( |2.7| , |2.8| ), will yield the bulk fields: a functional dependence for the phase-field and a 
differential equation for the stream function. 

Eq. ( p.8|) at 0(e°) and 0(e) reads respectively 

O(e ) : / o (0) = f(6 ) = 0^6 = 0,±1= const (3.5) 

0(e) : h{9) = -20x = =► e x = 0, (3.6) 

and iterating we get 

Q i = Vz > (3.7) 

Due to Eqs. (|3.5|) and (|3.7| ), 9 = ±1 to all orders, and, therefore, the (1 — 6 2 ) term in Eq. 
(277) does not enter this outer limit, whereas the viscosity contrast term in that equation 



becomes ±cV ip, depending on the phase. Hence, Eq. ( |2.7|) reads in this outer region 



e^ = (l±c)V 2 ^, (3.8) 



V 2 ^ = 0, ^ = (l±c)VV V*>0 (3.9) 



which implies 



except for c = 1. Note that we have recovered the sharp-interface Eq. (|2.1| ) in the e — » 
limit. For c = 1, Eq. ( |2.1| ) is still recovered in the +1 phase ( viscous fluid), whereas in the 
-1 phase (inviscid fluid) the stream function turns out to be constant in time to all orders. 
Although the inviscid fluid does not enter the problem in this limit (see Eq. |2.2| ), it still has 
a non-trivial dynamics, since the stream function in it must evolve to keep satisfying Eq. 
( |2.3p , and therefore, strictly speaking, we do not really get the right sharp-interface limit 
for c exactly equal to one. However, the model can be applied to physical high viscosity 
contrast pairs of fluids. We shall come back to this point in section |IV|. 
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B. Inner equations 



In turn, the interface boundary conditions for the stream function are given by the 
leading-order outer quantities ^^(O*) and — V'o^O - ). According to the matching 

conditions Eqs. (fyf) and (|3~4l ), these equal the inner ones \l/o,s(±oc>) and \E , i iP (+oo) — 
^i, p (— oo) respectively. Because of the specific structure of our phase-field equations, Eqs. 
( [2.7| , [2-8|) , we will need the first two orders in the inner version of Eq. ( p.8|) and the lowest 
one in that of Eq. ( |2.7| ) to get \]/o,s(=too), and the two first in Eq. ( |2.7p and the lowest in 
Eq. (|2.8|), to get \E f i i p(+oo) — ^ / i, p (— oo). Therefore, we compute the two first orders in both 
Eqs. ( |2.8| ) and (^.7[), by substituting the inner expansion Eq. ( |3T2"| ) in the inner (rescaled) 
equations (all whose terms are derived in the appendix): 

Eq. (pD up to 0(e) reads (see Eqs. (^, |A1^|ATT| , |AT9|) ) 



- ev n Q 0tP = f(Q ) + eQJ'iQo) + OiPP + e0 liPP + e(e , p ^ , s - ©o,^o, P )- 
Its O(e ) part, 

/(©o) + ©o,pp = 0, 

which, together with the boundary conditions specified by the matching (Eqs. 
with the outer expansion Eq. (|3.5|), gives the so-called kink solution: 



(3.10) 
(3.11) 

HID) 



©o = tanh -j= 
v2 



1 i 2 P 1/1 



Hence we find the 0o, s term to vanish, and Eq. (|3.10|) reads at 0(e): 

- v n Q , P = ©i/'(©o) + e liPP + ©o, P ^o, s 
As for Eq. Q, it reads, up to £>(±) (see Eqs. (^,|A^, [ATo1 , [A11 , |AT7| ) ) , 



(3.12) 



(3.13) 



it ^o,p P + -(*i,p P - ^o,p) + c{-(6 ^o,p)p + -[(©o*i, P )p + (©i*o,p)p - «©o*o,p]} 



+ -- 



1 1 



= 7 (i-eS) = o. 



From its O(^) part we know that 

^o,p(l + c @o) = const. 



(3.14) 



(3.15) 



Since \l/o,p has no correspondence with the outer expansion, it must vanish at infinity (Eq. 
( |3.4|) ). Then, we know the constant to be zero. Now, since the term in brackets vanishes 
only for c = 1, p — > — oo, we deduce that 



^o,p = 



We then put Eq. (|3~lp into Eq. ( gig ) at 0{\): 

^i,pp + c(e *i, p ) p = 



-|e , P 



(3.16) 



(3.17) 
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Finally, Eqs. ( |3.13|) and ( 3.17 ) will yield the macroscopic equations Eqs. (|2.3j ) and ( p.2| 
respectively: Eq. ( [3.13| ) can be rewritten in the form 



ie, = if(e ) + ^je, 



-e , P K + *o,.)- 



(3.18) 



We realize that LQ p equals the partial derivative with respect to p of Eq. (|3.11| ), which, 
in turn, vanishes. Hence, we write down the solvability condition 



/+oo 
(v n + ^ , s )Ol p dp = 
-oo 



(3.19) 



Using Eq. ( |3.16| ) we know that ^o,ps = = ^o,s P and can take ^ 0tS out of the integral as 
well as v n . Since the quantity left under the integral sign (0o,p) 2 is always positive, we find 
that v n + \l/o,s must vanish, and, matching with the outer expansion, we get Eq. ( |2.3j ) for 



On the other hand, integrating Eq. ([3.171) from p — > — oo to p — > +oo we get 

= -|©o - c6o^i, P + ai(s), 



(3.20) 



where ai(s) is an arbitrary function of s. Computing \E f l p (+oo) — ^ l p (— oo) and matching 
with the outer expansion Eq. ( |3.4| ) gives Eq. ( |2.2j ) for tp . This completes the sharp-interface 
limit. 



IV. FIRST ORDER CORRECTIONS TO THE SHARP-INTERFACE LIMIT 

In the phase field model the interface width and the convergence to the sharp interface 
limit is controlled by the small but finite value of the parameter e. Then, the question 
of which value of e is needed to accurately reproduce the actual Hele-Shaw dynamics for 
given values of the physical parameters B and c arises. This question can be qualitatively 
answered by noting the distinct roles played by e in the phase-field equations, Eqs. (|277] , p.8[) : 

The e factors appearing in e 2 V 2 6>, e 2 K(6)\V6\ and 7^7(0) (1 — # 2 ) ah stand for the 
interface thickness, and this is required to be small compared to the longitudinal length 
scale \k\~ l of the interface: e\k\ « 1. 

In contrast, the e in has nothing to do with the interface thickness (and we will 
therefore denote it by e from now on), but its aim is to ensure that the stream function is 
laplacian and not diffusive in the e — > limit, which commutes with the e —>■ one (the 
reader can convince himself of this by going through the limit again but now considering 
e of O(e )): e sets the time scale of the diffusion of the stream function through a given 
characteristic length of wavenumber k, nsp? ( see Eq- |5T8| ) , which must be much smaller 
than the characteristic growth rate of the interface M -1 , so that the stream function is 
slaved to the interface: 4^ << 1 ± c. We also realize that the viscosity contrast c can be 
arbitrarily raised, as long as e is correspondingly lowered. So our model is valid even for 
c — > 1, as long as this limit is taken formally after the e — > one. 

The e 2 in e 2 || represents the relaxation time of the phase field towards the steady kink 
solution (see Eq. |2.8| ), which must be kept well below the interface growth time M -1 for 
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the phase-field to remain close to the kink profile during the interface evolution: e 2 \uj\ << 1. 
This factor must be the same that the one in e 2 z- (Vtp x V#) in order to get the macroscopic 
equation Eq. (|2.3|) . In fact there are at least two distinct powers of e for this relaxation time 
(e and e 2 ) for which the right sharp-interface limit is achieved, and the corrections which we 
will compute would also be the same. 

To sum up, there are at least two independent small parameters (e and e) controlling the 
limit. When trying to approach macroscopic solutions by means of numerical integration of 
the phase-field equations, it is very convenient to vary them independently in order to save 
computing time, since both affect it [|23|| . 

A more quantitative answer to the question of the necessary values of e, e to get a given 
precision can be given by extending the asymptotic analysis of the previous section to first 
order in the interface thickness e considering e of 0(e). Thus, we will obtain a thin-interface 
model containing the corrections to the limit up to that order in e and e. 

According to the matching conditions Eqs. ( |3.3| ), the corrections to the interface bound- 
ary conditions for the stream function at first order in e, ^i iS (0 ± ) and ^i jr (0 + ) — ipi^(0~), 
are to be identified as terms in the expansion of \l/i iS (±oo) and ^2, P (+oo) — ^2, P (— oo) re- 
spectively. Now we will need the second order in Eq. ( p.8|) and the first in Eq. Q2.71) 
to compute 1 J r i )S (±oo), and the second in Eq. fl2.7p and the first in Eq. (|2.8p, to get 



^2,p(+oo) — ^2, P {— oo). Therefore, we must compute the next order both in Eqs. (|2.8|) and 
( f2.7[ ), but, first, we can still extract some information from the lower orders. 

On the one hand, we found that ^ 0jS = —v n . We put this into Eq. ( |3.13| ) to get the 
differential equation for 0i: 

9i/'(e ) + ei, pp = o (4.1) 



with boundary conditions coming from the matching Eq. ( |3.3|) with Eq. (|3.7|) @i(±oo) 
Oi iP (±oo) = and solution Oi = 0. 

The integral with respect to p of Eq. ( |3.20| ) is 



According to the matching Eq. ( |3.3|) , the p — > ±oo asymptotics of x 3/i(p) should consist 
of a finite term, ipi(0 ± ), and a diverging one, /?0o,r(O ± ). For vanishing viscosity contrast 
the last integral in Eq. Q4.2|) is pai(s) and clearly does not contribute to the finite term 
ipi (0 ). Then, since 0o is an odd function of p, its integral with respect to p will be even, 
and ipi(0 + )=ipi(0~), i.e., the fluid velocity normal to the interface will be continous on it. 
For non-zero values of c, however, one must compute the integrals in Eq. fl4.2p , find their 
p — > ±oo asymptotics, and identify ipii^) an d V ; o,r(0 =l= ). Requiring this latter quantity to 
be consistent with Eq. (|2.2|) for ipo, one fixes a±(s), and putting this back into the identified 
V'i(0 ± ) value, one finds 

^i(0 ± ) = -^{7 + c[^(0 + ) + ^o,r(0-)]}ln ^ + a 2 (s) (4.3) 

where 02(5) is another arbitrary function of s. This will give rise to a discontinuity in the 
fluid velocity: 
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<Mo + )-<MO") 



"^H + c[^o,r,(0 + ) + T/WCT)]} In i±£ 
-cv^K + c#o,r*(0 + ) + V"b,«(0-)]} + 0(c 3 ) 



(4.4) 



In order to fix d s a 2 (s), we compute the next order ((9(e 2 )) of Eq. (|2.8|) to get (see Eqs. 
(P HXT0| , |MT| , |XT9|) ): 



e 2 /'(e ) + e 2 , pp - p^e , P + *i,.e ,p = o (4.5) 

This has the same structure than Eq. ( p,13|) and an analogue solvability condition applies: 

/+oo 
*i,.©L^ = (4.6) 
-oo 



Substitution of the expression for ^ obtained by performing the integrals in Eq. (|4.2|) into 
this condition and subsequent computation of the resulting integral fixes d s a 2 (s) so that 



i~ [ y +c 2 ][1 



-7 + ±1 + ^ 



hi 



1+C, 



7s , ^0,rs(0 + ) + ^0,rs(0 ) lr 5__ , 2 2 3 



-][g=Fc+-^ + 0(<; 



Finally, to get ^ r 2 ,p(+ 00 ) — ^2,p( — oo) we need Eq 

^2,pp - K*l,p - <9 S V„ + C[(e ^2,p)p - K*l,p0O - d s V n Q ] + 3BpKK s Q 0tP = 

Integrating this from p ^ — oo to p ^ +00 we obtain 



2 ' 1 

(4.7) 

) at (£>(e )) (see Eqs. 

(4.8) 



[*2, 



+00 
pJ-oo 



/+00 
(1 + ce )*i,pdp - ^nH^ + c[e * 2l p]i~ = 
-00 



(4.9) 



where we have omitted integrals of odd functions of p. We use Eq. (|3.20| ) to rewrite the 
integrand of the remaining integral as — ^9 + ai(s). 6 is an odd function of p and does 
not contribute to the integral, whereas ai(s) gives rise to a divergent term of the type [p]t^. 
According to the matching Eq. (ft.4|), V ; i,r(0 ± ) corresponds to the finite part of ^2,p(±oo), 
so that we find 



W(o + ) - V»i,r(o-) = -c[^i, r (o + ) + <M(T)] 



(4.10) 



which will leave the jump of the normal derivative of the stream function across the interface 
unaffected at first order in the kink width. 

Putting Eqs. (glj,Iland HI for ^n JOLFlCl ) together, we get an effective sharp-interface 
model for the dynamics of the 6 = level-set up to first order in e and e: 



^ = (l±c)V 2 ^ 



Vv(o + )-Vv((T) 
^(o ± ) 



-r 



V2T S 



(4.11) 
(4.12) 

(4.13) 
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where T = 7 + c[ip r (0 + ) + ip r (0 )] is the weight of the vorticity defined in Eq. (|2~6|) evaluated 

up to 0(e) and <? ± (c) = 1 - 4, + (±1 + ^i) In f±£. 

Note that the desired corrections to the limiting equations Eqs. ( |2.1[ - |2.3|) in Eqs. ( 4.11 ) 
and ( |4.13| ) go as e and e respectively, and the fact that Eq. ( |2.2|) remains unaffected. Note as 
well that the correction in e appearing in Eq. (|4.13|) has nothing to do with an Allen-Cahn 
law. So the € 2 k(9)\V9\ term has cancelled this out even in the first order corrections. 



V. LINEAR DISPERSION RELATION UP TO FIRST ORDER IN THE 

INTERFACE THICKNESS 

In order to see how such corrections affect some relevant specific situation we compute 
the linear dispersion relation of a perturbation to the planar interface y(x) = Ae ut+lkx for 
Eqs. ( p~Hl - P~T3D . We make the ansatz 

i>{x,y) = a±A& ,t + ikx - q ^ (5.1) 

inspired by the actual Hele-Shaw result, where now the coefficient a± allows for distinct 
amplitudes in each phase for the stream function to satisfy the discontinuity in the normal 
velocities of Eq. fl4.13|) , whereas the decay length q± in the y-direction is set not only by the 
wavelength of the perturbation 2ir/k, but also by the diffusion length in Eq. fl4.11| ), which 
is also different in each phase. Thus, Eq. ( |4.11|) yields 



5± = l% ± , P ± = +^/l + p^ (5.2) 

In turn, taking into account that v n = ujAe wt+tkx and 7 = 2iAsign.(k)ujQe wt+tkx — where 
ojq = \k\(l — Bk 2 ) is the actual Hele-Shaw growth rate — , Eq. (|4.13|) fixes a± to be 

a ± = ^[l + e\k\V2g ± {c) P -±^-} (5.3) 
Finally, Eq. ( |4 . X 2| ) requires that the following dispersion relation is satisfied: 

- = n+ , % , [1 - e |fc| ^- +P+ ^ +a ~ cK + f- (1 + C)p - ] + 0(e 2 ) (5.4) 

(l+c)p- + (l-c)p+ [ 1 1 2 (1 - c)p+ + (1 + c)p_ J K ' K 1 

' e\k\V2^) + 0(c 2 ) + 0(e 2 ) (5.5) 



1 + f 6 



This consists of the well known Hele-Shaw growth rate multiplied by a factor smaller than 
1 carrying the corrections in e, e. We identify the conditions on e, e heuristically derived 
at the beginning of section [TV] to control how close this factor is to 1 and in general how 
close the stream function is to the actual Hele-Shaw one: euj/k 2 « 1 ± c (within p±) and 
e\k\ « 1 in Eqs. (I5.2H5.4T ); and the simplified version up to 0(c) eu/k 2 « 1 and e\k\ « 1 
in Eq. ( |5.5| ). The amplitude factor Eq. ( |5.3| ) can also be expanded in powers of c making 
use of Eq. (|5.5|) to find 
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a ± = -^( . . =F ce\k\V2) + 0(c 2 ) + 0(e 2 ). (5.6) 

Since these corrections have a stabilizing effect, they could affect the selection of the 
steady finger width. As a matter of fact, Ben Amar already showed that the d s (y ■ s) term 
of T s in Eq. ( |4.13| ) on its own was capable of selecting a finger width greater than 1/2 |2"2|| . 
Then, for small enough values of the physical surface tension (i.e., the physical selection 
mechanism), for which a width very close to 1/2 should be expected, this term could turn 
out to control the selection itself, so that an unexpected greater width could be obtained. 
Of course, this will not be the case if a sufficiently small value of the interface thickness e 
is used, so that the condition e\k\ << 1 is satisfied for the length scale set by the surface 
tension: e << \/B. 



VI. CONCLUSIONS 



We have introduced a phase-field model for Hele-Shaw flows with arbitrary viscosity 
contrast and shown it to yield the proper sharp- interface limit. We have actually found two 
independent small parameters (e and e) and three distinct conditions on them to control the 
convergence to the sharp-interface limit e, e — > 0. In particular, e must be lowered when c 
is increased. A thin-interface model, i.e. an effective sharp interface model keeping finite-e 
and -e effects, has been derived for the dynamics of the phase-field model up to first order 
in both of these parameters. This thin-interface model has then been used to explicitely 
compute the finite-e and -e corrections to the Hele-Shaw result for a specific situation such 
as the linear regime, thus suggesting that the single-finger width selection could also be 
affected by these finite-thickness effects. 



In the following paper |23[ we perform numerical simulations of the phase field model 



Eqs. (577,2.8), and we explicitly vary the two small parameters e and e independently. In this 



way we both control the simulation accuracy through the conditions mentioned to show how 
to reproduce the Hele-Shaw dynamics within this method, and explicitly check convergence 
in the interface thickness. 
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APPENDIX: 

Our goal here is to rescale the differential operators appearing in the phase field model 



Eqs. (|2.7|J2.8D . The first step will be to rewrite them in terms of the local coordinates defined 
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on the interface r and s. To do this, one must precisely define the curvilinear coordinate 
system and compute its so-called scale factors: 

Consider the 9 = level-set and its intrinsic coordinates s (arclength along it) and r 
(signed distance to it, positive for a point with 9 > 0), so that s x f = x x y. Let a be 
the angle going from x to s. Then k = a s is the 9 = level-set curvature. We introduce 
X, Y as the values of x, y for a point on the 9 = curve with a given value of s. By moving 
this point infinitessimaly along s we find that these values have changed in dY = (is sin a, 
dX = ds cos a. Consider also the coordinates of a point x, y with 9 ^ in terms of the values 
X, Y of its closest neighbour on the 9 = level-set and the signed distance between them. 
Taking into account that a is also the angle going from y to f, one finds x = X — rsina, 
y = Y + rcosa. Now one can compute the (positive defined) scale factors 

h 2 r = x 2 r + y 2 = l=^h r = l (Al) 
hi = x 2 s + y 2 = (X s — ra s cos a) 2 + (Y s — ra s sin a) 2 = 

= (cos a — rn cos a) 2 + (sin a — rn sin a) 2 = (1 — tk) 2 =^> h s = |1 — tk\ = 1 — rK (A2) 



Note that the last equality in Eq. (A"2) requires that rn < 1. In the inner region, where 
we make use of such formulae, this will hold as long as the interface thickness e is much 
smaller than the curvature radius at any point of the interface, i.e. not too far from the 
sharp- interface limit. Otherwise the present analysis would break down, because one could 
always find a point such that rn = 1, where h s would vanish, reflecting the fact that the 
change of coordinates has become ambigous in s. 

Then, the scale factors are used to express the differential operators in terms of r, s: 

Va = a r f H — s (A3) 

1 — rn 

V • a = (a r ) r + ~ n0j + (QS)s (a = a r f + a s s) (A4) 

1 — TK 

V 2 a = a rr - + + ™' a ' (A5) 

1 — TK {l—rK) z (1 — TK) 6 

Finally, one sets r = ep and expands in powers of e: 

(l-rn)- 1 = l + e P K + 0(e 2 ) (A6) 

One gets: 

Va = -^Apr + sA s [1 + epK + C(e 2 )] ( A7) 

V • a = \A r ) P + \{A S ) S - kA t ][1 + e P K + 0{t 2 )} (A8) 

V 2 a = ^A pp - - P k 2 A p + A ss + 0(e) (A9) 
This completes the rescaling. Capital letters denote fields written in the rescaled coordinates 



of the inner region. Any other quantity appearing in Eqs. ( |2.7| , p.8|) is derived from these 
ones. For instance we get: 
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da _ da , 
dt = di 



|x,y=const 



da _ v n , da . , 

— - v • Va = -— A p + — - v t A s + O e 
at e at 



(A10) 



where the partial (total) time derivative is computed keeping x, y (r, s) fixed, and -y is the 
velocity of the r, s frame respect to the x,y one, i.e., the interface velocity. Moreover 

z-(V^x V0) = ^{^e p [l + e P K + 0(e 2 )] - tf p 6 a [l + O(e)]} (All) 

v ■ (evil) = ^(m P ) P + (m s ) s - J«e* p (i + epK ) + o(e) = 

= ^m p ) P - - e KG* p - P K 2 m p + (6M> S ) S + 0(e)*„] + 0(e 2 ) (A12) 

The only terms left in Eqs. fl2.7 ,2.8) to compute are those containing 7(6?) and k(9). To 
construct them we will need the following quantities: 



Iea + ^[i + o(e)] 2 = +^ 



\ 



l + e 2 §§[l + C(e)] 2 
p 



B f 2 B 2 B f B 2 



(A13) 



(Note that B p > 0, since B is monotonic in p and we defined r to be positive for the 9 > 
phase) 



f (0) = W - H 



^ x Zi = ^ +Je ^ [1 + 0(e)] 

t X |V#l 1 + 6 



2 ie| 

2 6 2 „ 



+ C(e 3 



e 2 6 2 „ r B, 



(A14) 



We have termed this f(9) because it is indeed the unit vector normal to the 9 =const level- 
set on which it is computed. We similarly define s(9) = f(9) X z and k{9) = —V • 

generalized to construct = Bs(9) ■ V/c(0) + y ■ s(0): 



«(0) = V • 



V9_ 

W\ 



-« + e[( 



B, 



1 B 2 



2 v e 2 



(A15) 



Vk{9) 



'6, 



2 v e 2 - 



+ + 4-(^)ss + ^(|fU + 2p^ s ] + 0(e 2 )}[l + epn + 0(e 2 )] 



B, 



r\-^ + {—) sp -\{^) pp + 0{t)} 
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+ s{k 8 + e[3 P KK s - (|i) ss + )„] + 0(e 2 )} (A16) 
^ = fi{ Ks +e[3 P KK s - + i(|f )„] 

+y " S - e|jy " r + 0{e 2 ) = 1 + 0(e) (A17) 

We should still compute the product k(9)\V9\ appearing in Eq. (|2.8|), but instead we prefer 
to compute straightahead the sum V 2 8 + «(0)|V0| = V 2 0- V 2 6 + M ■ V| V0| = r(0) ■ V| V0|: 

V| V0| = f + i(S) p + 0( e) ] + s[^l + 0( e °)] (A18) 

v 2 + «(0)|w| = + _ ||e pp + 2|^e ps ] + 0(e) 

= % + + ^) - ^ + (A19) 



6 w p w p w p 
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